Ginzburg— Landau equations with consistent Langevin terms for 

nonuniform wires 



Jorge Berger 
Physics Department, Ort Bmude Colleqe, 
P. 0. Box 78, 21982 Karmiel, Israel t 

Abstract 

Many analyses based on the time-dependent Ginzburg-Landau model are not consistent with 
statistical mechanics, because thermal fluctuations are not taken correctly into account. We use the 
fluctuation-dissipation theorem in order to establish the appropriate size of the Langevin terms, and 
thus ensure the required consistency. Fluctuations of the electromagnetic potential are essential, 
even when we evaluate quantities that do not depend directly on it. Our method can be cast 
in gauge-invariant form. We perform numerous tests, and all the results are in agreement with 
statistical mechanics. We apply our method to evaluate paraconductivity of a superconducting 
wire. The Aslamazov-Larkin result is recovered as a limiting situation. Our method is numerically 
stable and the nonlinear term is easily included. We attempt a comparison between our numerical 
results and the available experimental data. Within an appropriate range of currents, phase slips 
occur, but we found no evidence for thermally activated phase slips. We studied the behavior of a 
moderate constriction. A constriction pins and enhances the occurrence of phase slips. 
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I. INTRODUCTION 

Almost a century ago Langevin [l| proposed a "complementary" stochastic force in order 
to enable the description of an individual Brownian particle in terms of Newton's second 
law. This approach gave origin to the field of stochastic differential equations. 

Among the many widely used stochastic differential equations, we will focus on the time- 
dependent Ginzburg-Landau equations (TDGL), which are used to describe a superconduc- 
tor out of equilibrium j^. The TDGL equations may be written in the form 

^hdttP = -[a + + |«|e'(^V - 27rA/$o) V , (1) 

adtA = {2hce/m)Re[ip*{iV - 27rA/%)ip] - {c^/A7r)V x V x A , (2) 

where i/j is the order parameter, A the electromagnetic vector potential, t the time, ^ the 
coherence length, $o = hc/2e the quantum of magnetic flux, a, (3 and 7 are material 
parameters, a is the normal electrical conductivity, c is the speed of light, e the absolute 
value of the electron charge, m the mass of a Cooper pair and the asterisk denotes complex 
conjugation. (3 and 7 are positive, whereas a is positive above the critical temperature and 
negative below it. The product |q;|^^ = li^ /2m does not depend on the temperature. In 
these equations the gauge choice eliminates the scalar potential. 

In order to cope with fluctuation problems, Schmid added a Langevin term to these 
equations [sl. Precisely, he added to the right-hand side of Eq. ([1]) a random function of 
position and time, f{r,t), such that 

(/*(r, t)/(r', t')) = 2m^kBT5{v - v')5{t - t') , (3) 

where (■ ■ ■ ) denotes ensemble average, VL is the volume of the sample, fc^ is Boltzmann's 
constant and T is the temperature. Since then the use of a Langevin term has become almost 
a standard practice in the numeric solution of the TDGL equations; some miscellaneous 
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10|. Generalizations of Eq. ([3]) to various classes of models 



examples are |4 
were reviewed in 

We shall be particularly interested in wires with nonuniform cross section. This interest 
is motivated by previous works [l2, [l^, which found that nonuniform superconducting 
wires may behave qualitatively differently than a uniform wire, especially near the transi- 
tion temperature. We are therefore interested in a numeric scheme that can enable us to 



investigate nonuniform wires in the region in which thermal fluctuations are important. For 
this purpose we would like to integrate the TDGL equations by means of Euler-Maruyama 
iterations. As we will discuss after Eq. (1161) . the existing literature does not provide an 
explicit procedure for evaluating the size of the fluctuations of the order parameter in each 
computational site; in particular, it is not mentioned that this size depends on the volumes 
of the computational cells. 



Moreover, addition of a Langevin term to Eq. ([2]) |10l . 1 151 . |16| , which can take into account 
fluctuations of the electromagnetic fleld, is the exception rather than the rule. We shall see 
that this omission is not self-consistent, even if only the values of the order parameter are of 
our concern. We are unaware of the initial motivation for including fluctuations in the order 



parameter only; a belated argument could be the result [17| that charges of an isotropic 
superfluid are irrelevant for ^ 1 K cm, but this result does not refer to the ID situation 
that we consider. 

An early review on fluctuations near superconducting phase transitions was written by 
Skocpol and Tinkham 181]. The voltage-current characteristic in quasi-one- dimensional su- 
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261 ] and reviewed |27l. l28] 



perconductors was extensively studied [19, 
several decades ago. The most salient feature is the presence of voltage steps, which are 
atributed to "phase slips" in which the order parameter vanishes at a point and releases 
a winding from its phase. In recent years there has been a renewed interest in this topic 
29, 3^, 31, [32, 33I, 34], fuelled by the availability of long nanoscopic superconducting wires. 
Among other features, it has been possible to investigate the existence of quantum phase 
slips and the behavior under flxed voltage. 

In this article we build the equations for the evolution of a nonuniform superconducting 
wire, as appropriate for a flnite difference numeric treatment. For this purpose we flrst 
discretize the TDGL thermodynamic potential and then invoke the fluctuation-dissipation 
theorem 351] to evaluate the self-consistent Langevin terms. 

Our method is developed in Sec. [Tll Expressions for the noise term at each computational 
element are obtained and we flnd a convenient way to determine how the Johnson noise 
affects the order parameter. Section IIIII compares the results of our method with those 
of statistical mechanics for simple situations in which analytic results can be obtained. 
In Sec. IIVI our method is used for situations in which analytic results are not available. 
Appendix El provides numerical details that contributed to the success of our algorithm. 
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II. OUR METHOD 



A. The fluctuation-dissipation theorem 

For our purposes, the theorem may be stated as follows. We consider a system with 
energy G that depends on several variables, one of which is x. Let the system be in thermal 
contact with a heat bath at temperature T and let x be a classical variable, i.e., a real 
function of time. Let x obey on the average the macroscopic equation 

dx/dt = -TdG/dx . (4) 

Let r be a macroscopically short period of time. Then the microscopic change of x, x{t + 
t) — x(t), will be —TrdG/dx on the average and its variance will be 

= 2TkBTT , (5) 

where 77 = x{t + r) — x{t) + TrdG/dx. 

The physical idea behind this theorem is the same as in the original Langevin article, i.e., 
the same interactions of the variable x with the heat bath, that give origin to its microscopic 
fluctuations, are those that lead to its macroscopic evolution in Eq. (jl]). 

Several assumptions are necessary for the validity of this theorem. One of them is that the 
distribution of the coordinates of the heat bath that interact with x remain described by the 
canonical ensemble. There are also conditions on the size of r; it has to be macroscopically 
small, quantitatively expressed by r ^ ksT /V{dG / dxY and r ^ 1/Vd'^G / dx"^ . On the 
other hand, r has to be microscopically large, explicitly, much larger than the autocorrelation 
time of the Langevin term and also large enough to allow for a classical treatment of x. From 
the numerical point of view, taking a value of r that is too small may be a waste of resources 
but, in principle, should not affect the results. Indeed, if we take a time step t/K rather 
than r, the average and the variance of the change in x will decrease by the factor i^, so 
that, if 7] has Gaussian distribution, after K steps the ensemble of results will be the same 
as if the time step were r. 
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B. Ginzburg-Landau energy 



In this model the energy of a superconducting wire can be written as 

G = j[a\ij\^ + P\tP\y2+\a\^^\{T^ -2nA/<^o)ij\^ + IA/cw]dV , (6) 

where the integral is over the volume of the wire. Here / is the current that flows along 
the wire, w is its cross section, and the last term in the integrand takes into account the 
energy provided by the power source. We consider a wire with a cross section that is much 
smaller than the square of the magnetic penetration depth, so that the induced magnetic 
field can be neglected. In this situation we are also allowed to take the electromagnetic 
potential A in the direction of the wire. Likewise, we will assume that the width of the 
wire is small compared to the coherence length and varies slowly with the arc length; for 
the superconducting-insulator boundary condition these assumptions imply that V may be 
replaced with the derivative with respect to the arc length. 

Let us now discretize this energy. We denote the length of the wire by L and divide it 
into segments of length L/N. Each segment has to be sufficiently short, so that the order 
parameter, the elecromagnetic potential and the cross section may be considered to remain 
uniform within it, and we denote their values by ipk, and Wk for the k}'^ segment. Defining 
the grid-dependent dimensionless variables ^ = NC,/L and A = 2tiLA/N^q, we may write 



N 



2 



~e [1(1+ iAk)'iPk - 'ipk-i? + i(-i + iAu)'ipk + v-fc+ii'] } + ^ SZ^'^ • ('^) 



The intuitive meaning of is "the number of consecutive segments over which the order 
parameter remains approximately uniform." In numeric calculations ^ has to be of the order 
of 1; if is too small the discretized model does not represent the physical sample and, if it 
is too big, the calculation becomes inefficient. 

Discretization of the local terms in Eq. ([6]) is quite natural, but discretization of the term 
with the gradient is somewhat arbitrary. Since the term (iV — 27rA/$o)V' requires the values 
of the vector potential and of the derivative of ip at the same site, it is usually convenient 
to use a staggered approach, in which A and if) are defined at alternating overlapping sites. 
However, since in Eq. |(iV - 27rA/$o)V'P is muhiplied by the cross section of the wire, 
which may not be uniform, the use of overlapping sites is problematic. It is tempting to 



replace in Eq. ([7]) the term in square brackets with l^/^fc+i— V^fc-i + 2«Afc'?/'fcp/2. Unfortunately, 
even in the limit ^ 1, this form would lead to two separate sublattices {k odd and k 
even) such that there is no penalty for changing their relative phases; this feature does not 
represent the original continuous model. 

We decompose ipk into real variables ipk = Uk + ivk- The derivative of the discretized 
energy with respect to Uk may be written as 

dCjf Lxv 

— = -^Re |2a^fc + 2/5|7/'fc| Vfc + 1"!^' [2(1 + i^^fc - (1 + - (1 - 

w^^[wk-i{tpk - {I - iAk-i)'ipk-i) + Wk+iitpk - {'i- + iAk+i)tpk+i)] \- (8) 



The derivative with respect to Vk involves the imaginary part of the same expression, which 
is actually 2dG/dijjl. 

The derivative of the energy with respect to Ak is 

—~- = n ^ ? - Vk-^l + 2AkUk)Uk - {Uk+l - Uk-l - 2AkVk)Vk\ . (9) 

dAk 2tcc L 

Noting that the term in square brackets is proportional to a discretized version of the 
superconducting current density and substituting |q;|^^ = !f/2m, Eq. ([9]) reduces to 

dG In^o 



dAk 27rc 
where 1^ is the normal current along the wire. 



(10) 



C. Evolution equations 

Comparison of Ohm's law with Eq. ffTOj) gives the macroscopic equation 

dAk/dt = -TA,kdG/dAk , (11) 

with 

= NW^k ■ ^^^^ 
Equation ([8]) has to be compared with the evolution of the order parameter in a wire 



14| . Using superconductor- insulator boundary conditions and our notation, this evolution 



is given by 

^hdtij = -{a + (3\^P\^ + \a\e[{tds - 2TiA/%f + i{w\s)/w{s)){id, - 27rA/$o)]}^ , (13) 



where s denotes the arc length. Assuming that w, and A can be differentiated twice, 
Eq. f[T3|) is obtained as the hmit of 



duk/dt = -T^^kdG/duk , (14) 
and an analogous equation for Vk, with 

= ■ 

The macroscopic evolution of Uk and Vk can be described by the complex equation dipk/dt = 
-2T^,kdG/drk- 

We close this section with a summary of our basic numerical algorithm. Let tpk and Ak 
be known at a certain instant. Then, after a time step r their new values will be 

dG 

di^k 
dG 

Ak ^ Ak- Tamt—y- + Va , (16) 
oAk 

where t]u, rjy and r]A are random numbers with gaussian distribution, zero average and 
variances {riD = (rjl) = 2T^^kkBTT, {ri\) = 2V A^kkBTr, with the derivatives of G given by 
Eqs. (EJ-Q and the F's given by Eqs. (fT2l) and ffTSj) . It should be emphasized that TA,k 
and depend on the length and the cross section of the k}'^ segment. 

In order to use algorithm ( fT6l) . (tiI^^^a) must be known. Nevertheless, practically all the 
literature limits itself to mentioning the continuum-style relation ([3]). Neither the expressions 
nor the procedure for evaluating {rjl „ a) ^-re provided. To our knowledge, the only exceptions 
are 4| and The expressions provided in these references neither coincide with those we 
have obtained nor among themselves. According to us, the expression in Eq. (16) of ^ 
should vanish for a truUy infinitely long sample; expression (18) in seems to depend on 
;he volume of the entire sample rather than the volume of the computational cell. Reference 
Gj just quotes Q and 'although the size of the noise is essentially regarded as a free 
computational parameter), and jol relies on 6]. Many studies (e.g. [s]) do not require the 
variance of riu,v,Ai because algorithm f|T6|) is actually not used. A possibility for avoiding use 
of (fT6|) is working in Fourier space, but this approach is not useful for a nonuniform wire. 
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D. A gauge- invariant method 



For a one dimensional wire the electromagnetic potential can be gauged out of the evo- 
lution equation by means of an appropriate transformation. Defining a gauge function 

W(s) = expf^ rA{s')ds'] (17) 



and the gauge-invariant order parameter ipi^) = ^{s)il){s), the term |(iV — 27rA/$o)'V^| in 
Eq. ([6]) reduces to \d%lj/ds\. We discretize lA as 

/ fc-l/2 ^ \ 

= exp \ i Yl^^n ' (IS) 

where the sum up to a half-integer number is understood to contain the term A^, but 
multiplied by i. Once is defined, we can write ipk = WfeV'fc- A discretized version of the 
energy can then be written as 

i^l'K + ^fc+OlV-fc+i - + ^ > (19) 

where Wk is understood as for = and k = N + 1. 

Although ipk is physically more meaningful than ipk {\d'4)/ds\ is proportional to the ve- 
locity), it is not a "canonical" variable, i.e., dipk/dt ^ —2VdG / difjl- Instead, we have 

~ fc— 1/2 

Wfc— - + #fe2^— -. (20) 



Noting that 



dt dt ^ dt 



the first term in Eq. fl20l) takes the value —2T^^kdG/dip'l,, i.e., its contribution is what we 
would obtain if ipk were a canonical variable. In order to obtain the evolution of Aj, we 
express the energy in Eq. ( |T9l) in terms of ip rather than ip and then use Eq. (ITT]) . This gives 

^ = -r^j + ^"^^^ Im (wj_i + «;j)Vi;_i^j + (wj + Wj+OV'JV'i+i j , (22) 

which is the discrete version of Ohm's law required by Eq. ( |T9l) . Substituting these results 
into Eq. (!20!) we obtain the macroscopic evolution of ip. This evolution can be followed 
without any knowledge about A. 
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We investigate now the influence of fluctuations, tpk is the same as tpk, rotated in the 
complex plane; therefore, its fluctuations arise from those of ipk and from those in the angle 
of rotation. Since the size of the fluctuations of the real and the imaginary part of ipk are 
the same and since they are not correlated, their contribution to the fluctuations of ipk has 
the same distribution as rju + irj^. The influence of the fluctuations in the angle of rotation 
can be taken into account as follows. When integrating the evolution of ^pk in Eq- ( |20l) we 
require the macroscopic change of Aj, rdAj/dt. To this macroscopic change, we have to add 
the fluctuation of Aj during the period r, rjA, which is still described by {ri\) = 2V AjkBTr. 

Our algorithm for the evolution of ipk consisted of compound steps. The flrst stage of a 
step, that corresponds to the flrst term in Eq. ( !20l) . looks the same as the evolution of ipk 
in Eq. (fT6l) . with A^ set to zero. The second stage takes care of the last term in Eq. (l20ll ; 
instead of a standard Euler step we multiply by a phase, i.e.. 



where the meaning of the upper limit is that the argument of the exponential for j = k is 
divided by 2. To flrst order in r and rjA,j this method is equivalent to an Euler step, but, 
since it keeps unchanged, is numerically more stable. 

E. Dimensionless parameters 

Besides ^ and A, that were deflned before Eq. ([7]), it will be convenient to deflne the 
following dimensionless quantities, a = a/ksT provides a conveniently normalized value 
of the condensation energy. L = L{2mkBTy'^/h is the length of the wire divided by the 
"thermal wavelength," except for a factor ^/^^. We also deflne R = Ae'^'yL/hcrwQ, where wo 
is some typical cross section of the wire. For a wire of uniform cross section, R is the normal 
resistance of the wire multiplied by 2717 and divided by the quantum of resistance, /i/(2e)^. 
R is proportional to the ratio between the relaxation times of the order parameter and the 
electromagnetic potential, f = ar/'yh is the effective size of the iteration steps; it is not 
dictated by the physics of the problem and we expect that, provided that f is signiflcantly 
smaller than 1 and signiflcantly larger than the inverse of the number of iterations, the 
results should not depend on it. Some additional quantities, which are useful when /3 7^ 0, 




(23) 



will be defined in Sec. IIV C 1[ 
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III. TOY MODELS 



In order to test our method, we would like to compare its results against quantities that 
can be evaluated exactly. For this purpose, the following sections consider simplified models. 



A. Uniform A 

We drop the last term in Eq. ([7j), that causes the energy to be unbounded from below, 
and also set /3 = 0, so that the problem becomes linear in the order parameter. In order 
to keep G bounded from below, we must set a > 0. As additional simplifications we will 
consider a uniform wire Wk = Wq for k = 1, . . . , N, and a uniform value = A for all the 
segments. With these simplifications the energy becomes 

G = ^ E I l^'^l' + I + ^^)^'^ - ^'^-il' + + ^^)^'^ + } • (24) 

We also have to specify boundary conditions at the extremes of the wire. Let us consider 
periodic boundary conditions, ipQ = ip^, i^N+i = V'l- 

For a fixed value of A, G can be diagonalized by a Fourier transformation. We define 
ifn = N~'^ Yl V^fcC"^'^*'^"'/^ and, for simplicity of notation, we will take even values of N. The 
inverse transformation can be written as 

N/2 

i^k= E ^.e^-^'^"/^. (25) 

n=-N/2+l 



Substituting Eq. ([251) into Eq. ([24]) we obtain 

N/2 

G = aLwo J2 IVnl"^ {l + ^\A + sm2TTn/Nf + Asm* nn/N]^ . (26) 

n=-N/2+l 

The ensemble average of any quantity Q is given by 

(Q) = J QiWn}; i)e-^/'=^^dV-iv/2+i ■ ■ ■ c^Viv/2cii/ J e-^/^^^ciV-iv/2+1 ■ ■ ■ d^^N/2dA , 

(27) 

where A varies over R and cpn over R^. If Q is a polynomial, the integrals over ipn are easily 
evaluated. For example, the denominator in Eq. 0271) (the partition function) is 

{iiksT / aLwo)^ dA 



;i + e(i2 + 4)) nn5-W2(l + ^~'[(^ + sin 27rn/iV)2 + 4 sin^ vrn/iV]) 



(28) 
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This last integral may be either evaluated numerically or collecting the residues at A = 
sin 27m /N + + 4 sin^ 7m /Ny/^. 

Note that, since the dependence on A factors out in the integrand of Eq. (1281) . it con- 
tributes just a multiplicative constant to the partition function. Therefore, A is a degree 
of freedom that makes no contribution to the heat capacity. For the same reason, the 
expectation value for any function of A is independent of the temperature. 

In our calculations using algorithm f|T6|) we started from Uk = = A = and let these 
variables build up from fluctuations. We performed A^reiax iterations to achieve a "typical" 
situation and then averaged the quantities of interest during A^^av additional iterations. 

Figure [1] compares the averages of several quantities, as obtained with our method for 

= 4, with those expected from Eq. (!27|) . for a wide range of values of ^/L. If the 
fluctuations of A are ignored, we obtain the results in the lower part. As expected, ignoring 
the fluctuations of A yields a variance of A (not shown in the graph) that is much lower 
than its ensemble value, but we see in the lower part that also the averages of quantities 
that do not directly involve A turn out to be incorrect. 

The upper panel in Fig. [2] compares the averages (Iv^nP) obtained with our method with 
those expected from Eq. fl27|) for A^ = 10, and the lower panel repeats the comparison while 
ignoring fluctuations of A. We see again that the symbols in the lower panel deviate from 
Eq. ( l27l) . but this time the disagreement is smaller than in the case A^ = 4. The reason is 
that the variable A is driven by the other variables Ui and vf, when C,/L is sufficiently large, 
the variance of A is practically independent of whether or not its fluctuations are considered. 
This is a feature of our toy model, in which there are 2N degrees of freedom for the order 
parameter and only one for A; in the original model with A^ degrees of freedom for Ai the 
influence of the fluctuations of the electromagnetic fleld should not decrease with A^. 

B. Gauge-invariant formalism 

We still take I = f3 = and uniform cross section = Wq, but use of the formalism 
developed in Sec. IIIDI enables us to release the requirement of uniform A. We consider two 
boundary conditions. One of them is V'(O) = 'ipiL) = 0, which corresponds to contacts with a 
conductor that strongly suppresses superconductivity. In a discrete model, these conditions 
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FIG. 1: Averages of several quantities, as functions of the coherence length. The symbols corre- 
spond to values calculated using our method (upper part) and the lines are ensemble averages. 
The parameters used in the calculation were = 4, f = 3 x 10"^, R = 16000, A'av = 2 x 10^, 
A'rclax = 6 X lO"^. The parameter d factors out. • dLwo{\(po\'^) + 0.25 (the constant 0.25 was 
added for visibility); * dLi^odv^-iP) + 0.25 and dLwQ{\ipi\'^) + 0.25; + dLit;o(|v'2p) — 0.25; A 
{dLwo)'^{\ipo\'^); □ {dLwo)'^{\ipi\'^); {dLwof{\(po\'^\ipi\'^); x 0.5 logio(^^)- Upper part: fluctua- 
tions of A were included; lower part: fluctuations of A were ignored. We have drawn an oblique 
straight line that separates the "upper" and "lower" parts. 

are approximated by ipo = —ipi and ipN+i = —i'N- In this case we write 

7 , n(k -l/2)7r , , 

= 5^ sin ^ ^' ' (29) 

n=l 



with 



, 2-5nN^j . n{k -1/2)tx 
^n= ^ l^^ksm , (30) 

k=l 



where 6 is Kroneker's symbol. Expression ([T9l) then assumes the diagonal form 

G = aLwoJ2i^ + S^^)MU- + 2i'sm'^] . (31) 

n=l ^ ^ 

The other boundary condition we consider is that in which dip/ds vanishes, which cor- 
responds to contacts with an insulator. In a discrete model, we require ipQ = ipi and 
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FIG. 2: Variances of the Fourier components of the order parameter, as functions of the coherence 
length, for the case of 10 segments. The curve for n = 5 sits at its true position; for visibihty, the 
other curves have been raised by (5 — |n|)/10. Besides N, the parameters are the same as in Fig. [TJ 
• n = 0; ★ n = ±1; A n = ±2; □ n = ±3; n = ±4; x n = 5. Upper panel: fluctuations of A 
included; lower panel: fluctuations of A ignored. For visibility, only the region above 0.6 is shown. 



ipN+i = i^N- In this case we write 

i'.-Y.A^o.'-^^^^^ (32) 

n=0 

with 

^n=^^2^^kCOS (33) 

k=l 



and Eq. (fT9|) reduces to 



f /I 

I n=l ^ 



^-n^^ir (34) 



Figure [3] shows values of (Iv^^P) and (Iv^^P) obtained from Eqs. (1501) and after 
following the evolution of according to Sec. Ill D[ We see that, in spite of the fact that 
the vector potential has been gauged out and the macroscopic evolution of ip can be 
described in closed form by means of Eqs. (!20|) and (l22l) . the fluctuations of Ak still remain 
essential. 
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FIG. 3: Variances of the Fourier components (p^ and c/?^ of the gauge-invariant order parameter ■0, 
for = 6. The dashing Unes describe (Iv^nP) according to Eq. ([HT]) (-0 = at the boundaries), and 
have been raised by 0.4(6 — n). The sohd hnes describe (Iv'^P) according to Eq. (j34p (dip/ds = at 
the boundaries), and have been raised by 0.4(5.5 — n). Note that (ly'oP) does not depend on ^. The 
symbols describe the values obtained using the method developed in Sec. IIIDI • (Iv'.^P), n odd; * 
ilfnl'^), n even; A (Iv^^iP)! n odd; □ n even. Besides = 6 and f = 10~^, the parameters 

are the same as in Fig. [TJ Upper panel: fluctuations of A included; lower panel: fluctuations of A 
ignored. In order to avoid clutter, only the cases n = 0, n = 1 and n = 4 are shown in the lower 
panel. 

Inspection of the lower panel in Fig. [3] (including 7 additional curves that are not shown) 
suggests the following trends: (i) omission of the fluctuations of A yields values of (|v?nP) 
that are larger than the true values for small n and smaller than the true values for large 
n, i.e., the order parameter appears to be flatter than it should; (ii) the largest deviations 
of (Iv^nP) from its true values occur for coherence lengths of the order of the length of the 
wire; for higher harmonics, the largest deviation tends to occur at smaller coherence length. 

Trend (i) might be understood if we consider the mechanism by which fluctuations of 
Ak influence the order parameter. For a bulk superconductor, the currents generated by 
these fluctuations induce magnetic fields that are felt by the order parameter; this is the 
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mechanism considered, e.g., in 36|. For a very thin wire, the induced magnetic field is 
neghgible. However, since the total current has to be uniform along the wire, the Johnson 
current forces a superconducting counter-current, which in turn forces a phase gradient. 
Therefore, fiuctuations of force the order parameter to be rougher than it would be in 
their absence. Note that at a very small scale, which is beyond the scope of our model, 
electroneutrality is not a reahstic assumption. Trend (ii) may be qualitatively understood, 
since higher harmonics probe shorter lengths. 



C. The continuous limit 



We consider now the limit N oo. In this limit Ta 0, and we might expect that the 
fluctuations of Ak become unimportant. As in the previous section, we drop the nonlinear 
term, take a uniform cross section, and use the gauge- invariant order parameter ip. This time 
we will consider the periodic boundary condition ipN+i = "ipi, which physically corresponds 
to a ring that encloses a magnetic flux that equals an integer number of quanta. The case 
of a ring that encloses a non-integer flux will be considered elsewhere. For this boundary 
condition we cannot take / = 0; instead, / becomes a Lagrange multiplier. Numerically, the 
effect of this current amounts to the substraction of a uniform term from every Aj., so that 
^y4fc(t) remains constant. Since oc N, dip/dt diverges in the limit —>■ oo; in order to 
overcome this difficulty we used the method described in Appendix lA 1[ 

Defining = ^ V'/cC"^'^*^"'^^, the statistical average of (|<^nP) can be obtained by 
setting A — > in Eq. fl26|) . From here we obtain (|<^nP) = ksT/aLwo^l + AC,"^ sin^ {nn / N)) . 
Table [T] compares this statistical average with the values obtained with and without fiuctu- 
ations of Ak for several values of A^ and n. For the parameters considered, we see that the 
larger the number of elements into which we divide the sample, the more severe the error of 
the results obtained when the fiuctuations of Ak are ignored. 

Since the disparity between the rows a and c in Table [T] grows as the grid resolution is 
increased, one should ask whether this effect is due to the necessity of resolving small-scale 
features, or just due to the lager number of computational cells. In order to answer this 
question, we repeated our calculations for A^ = 32, but increased R, L/C, and by a factor 
of 4. In principle, one should recover the results obtained for A^ = 8, and this is indeed 
found when the fluctuations in Ak are taken into account. In row c, we recover the values 
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TABLE I: lO^aLwod'^nP) for several values of the number of elements N and the wave number 
n. We have taken a coherence length equal to the length of the sample, and the other parameters 
are as in Fig. [TJ 
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iV = 16 


iV = 32 
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10771 


47427 


10000 


11249 


71936 


10000 


9750 
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78 


49 


66 


66 


37 


64 


63 


34 
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46 
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29 
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28 


16 
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39 


39 


39 


19 


20 


14 


17 


17 


10 


16 


16 
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"Statistical average 

^Evaluated with the method of Sec. IIIICI 
'^Method of Sec. IIII CI without fluctuations in Ak 

for N = 8 for n ^ (i.e., these harmonics behave as "intensive" quantities), but the value 
obtained for = is roughly multiplied by 4 (behaves as "extensive"). 

It might be suspected that our numerical method breaks down for large A^. This is very 
unreasonable. First, practically the same results were obtained for different values of r. 
Second, if the results without fluctuations in Ak were numerically wrong, one would need a 
miracle to correct the results by taking these fluctuations into account. 

D. Nonuniform ring 

We finally consider a case with nonuniform cross section. We use the gauge-invariant 
description in Eq. flTI?]) . with j3 = and consider a ring with just A^ = 3 elements. For 
simplicity, we also take ^ = L/3. We take periodic boundary conditions, as in Sec. IIII C[ 
We denote the cross sections by Wi = ^Wq and W2 = = Wq. This is an artificial model 
and clearly we cannot regard its boundary as a smooth surface, but the important feature 
for the present purpose is that it is a self-consistent model. 

The energy can be diagonalized by a linear transformation ipi = (fi + (f2 + f3, V^2,3 = 
- (/i/2)((^2 + V3) ± ^iV2 - <^3), with u = i((4 + 12/i + 7fi^ + fi^)/{7 + This 
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FIG. 4: Statistical averages of several fluctuating quantities as functions of uniformity. For a 
ring with uniform cross section, /x = 1. The lines were obtained using statistical mechanics, the 
empty symbols were evaluated using fluctuations with the sizes dictated by Eqs. ([12]) and ([15]) 
and the fllled symbols were evaluated using fluctuations of the same size for the three elements. 
A aLwoilipil"^); □ aLii;o(|(^2,3p) i^i,2,3 as defined in Sec. [mD|) : o {G)/2kBT; o {G'^)^''^ /2kBT . 
Besides = 3 and = the parameters are the same as in Fig. [TJ 

transformation gives G = {aLwo/6){2 + fi){2\(pi\'^ + (2 + 5/i + yU^)(|v92p + Iv^sP))- 

As in Sec. IIII C[ the current becomes a Lagrange multiplier. Since ^ oc w^^, this time 
we have to substract a uniform term from every WkAk in order to keep a constant value of 

Figure m shows the statistical averages of |v?i,2,3p and of some powers of the energy that 
are obtained for this model using the method developed in Sec. [Til For comparison, we 
also evaluated these averages using fluctuation sizes that are independent of the local cross 
section; for the purpose of this comparison we replaced Wk in Eq. ([T^ and Eq. ([T^ with 
the geometric average of the cross section, fi^^^vjQ. 



E. Paraconductivity under constant electric field 

Tucker and Halperin 



2l\ evaluated the supercurrent along a thin uniform wire when 
a longitudinal electric field is applied. The strongest assumption of their model is that 
the electric field remains constant in time and position, and is not influenced by ip. An 
unphysical feature of this model is that the total current is not uniform along the wire. The 
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FIG. 5: Average of the superconducting current as a function of the apphed voltage under the 
assumption that the electric field is constant and uniform, in spite of the fluctuations of the 
superconducting order parameter. The line corresponds to Eq. (j35p . In all cases f = 3 x 10~^, 
I = 21 and, unless stated otherwise, A^'av = 2 x 10^, A'^rciax = 6 x 10'' and a = 10. A ^ = L, q = 4; 

□ ^ = 0.3L; O ^ = L; X ^ = 3L, iVav = 3 X 10^ iVrelax = 9 X 10^. 



supercurrent is defined by its spatial average, as in Eq. fIBip . An additional assumption 
was that the length of the wire is infinite. The main effort of Ref. 2l|] was focused on the 



influence of the term Plipl"^ in Eq. (Q, but at this stage we keep /? = 0. Their result was 



'sc 



eksT 



dy y^^'^ exp 



2/3 



y - 



12 



(35) 



where V is the applied voltage and Vq = aL/'je^. /sc is an increasing function of V; for 
V ^ Vq, /sc approaches 2ekBT / ^/Sh = LlbeksT/h. 

In order to reproduce this model, we follow just the evolution of ipk, whereas Ak{t) 
follows from dA^/dt = —2eV/Nh. We take periodic boundary conditions for ip, i.e. ipN+i = 
exp{—2ieVt / fijipi, and the condition of infinite length requires L ^ ^. The results are shown 
in Fig. [5] for several cases in the range 0.3 < < 3L. There is good agreement with the 
prediction of Ref. 2l| but, surprisingly, the requirement of large L seems to be unnecessary. 
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IV. APPLICATION PARACONDUCTIVITY 



A. Voltage current characteristic of a uniform wire above the critical temperature 

1. Fixed voltage 

We are now in a position that enables us to evaluate self-consistently the average current 
(J) along a superconducting wire when a fixed voltage V is applied. The procedure is similar 
to that of Sec. IIIIEl but this time the evolution of is given by Eq. f|T6l) . The instantaneous 
total current I has to be such that in every step ^ decreases by 2eVT/h. (/) is obtained 
as the average of / over N^v iteration steps. 

The dashing lines in Fig. [6] show the values obtained for {I){V) in the range 10^^ < 
Ih/eksT < 10^/^ and 40 < ^ < 25000. For ^ = 40 these values are scarcely below the 
currents that would be present in a normal material (in a logarithmic scale); this means that 
supercurrents are small compared to the normal currents. For larger values of R, Ohm's 
law is approached in the large current extreme, whereas in the low current extreme (/) is 
not proportional to V. It might be argued that (/) could become proportional to V for 
still lower currents; we did not attempt to evaluate (J) for lower currents, because Eq. (]B5P 
predicts fiuctuations of Isc^/eksT of the order of \f^jL for a single step and, since the 
number of independent values of Jsc is of the order of A^av^" ~ 10^, we should at best expect 
{I)h/ekBT to be accurate to the order of 10^^. 

2. Fixed current 

Another experimental situation is prescribed by keeping a fixed current J. In this case 
we are interested in the average value of the voltage. The evolution of ip and A is the same 
as in Sec. IIII CI except for an additional term RfhI/2NaekBT in the increment of each 
Ak in every step. The boundary conditions are dictated by the experimental setup; in the 
absence of experimental details we shall just assume that the coupling between ipi and ipo 
(and between ip^- and ipN+i) can be neglected. In order to avoid a significant infiuence of 
the boundary conditions at the "contacts," we adopted a "four-probe technique," i.e., the 
voltage was measured between the segments 1 + A^cont and — A^^cont and the result was 
multiplied by A^/(A^ - 2Ncont)- 
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-1 -0.5 0.5 1 1.5 

logio (I*/ekBT) 

FIG. 6: Voltage-current characteristic. Filled symbols (dashing lines) stand for {I){V) and empty 
symbols (continuous lines) for {V){I). The straight lines correspond to Ohm's law in the normal 
state. Unless stated otherwise, N = 21, A'^cont = 3, iVav = 2 x 10^, A'rciax = 6 x lO'^, d = 10, 
^ = L/\/To and f = 10~^. The values of R are marked next to the curves in the large-current 
regime. The lines with the rhombs correspond to i? = 1000 and a = 3; for visibility, they have been 
lowered by half a unit. For R = 25000, several points of {I){V) were evaluated using f = 3 x 10^^ 
and several points of {V){I) were evaluated using f = 3 x 10-6, iVav = 10^ iVreiax = 3 X 10^; thsse 
changes of computational parameters do not produce noticeable changes in this figure. 



The average voltage drop (V) is given by 



^ 5Z Aifc , (36) 



where AAk is the increment of Ak in a step and the sum is over the segments 1 + A^^cont < 
k < N — A^cont and also over the A"av iteration steps. 

The continuous lines in Fig. [H] show our results for {V){I). For large and intermediate 
currents {V){I) coincides with (/)(V^) in the scale of the graph, but for low currents we see 
that they are different curves. This difference reflects the fact that a situation in which / 
is kept fixed and V fluctuates is not equivalent to a situation in which V is kept fixed and 
/ fluctuates. This time we do reach a small current regime where (V) is proportional to /, 
with an effective resistance that is smaller than L/wqct. 

The difference between (V")(/) and {I){V) has been dramatically observed in recent ex- 



periments 



32 



33| 
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B. Low current limit 



Let us denote by a' the contribution of the supercurrent to the electric conductivity. 
Taking the hmit ^ ^ in Eq. fl5Sl) we obtain 

a' = aR^/8aL . (37) 



This resuh was initially obtained by Aslamazov and Larkin 



37| and the value of a' according 



to Eq. (1371) will be denoted by ctal- It can also be written as 

ctal = RaL\^/2Lf = Ra/Sa^^'^L = Ke-^^^ , (38) 

where e = (T — Tc)/Tc, is the critical temperature and K = ne'^^{0)/16hwo. The last 
form of Eq. ( l38l) neglects terms of order e~^/^ and assumes a = 8/cb(T — Tcj'y/ir. 

We will now evaluate a' taking the electromagnetic fluctuations into account. This will 
be done by fixing a small current I, evaluating the average voltage (V^) as discussed in the 
previous section, and then a' will be given by 

^ RhI/ekBT _ 
a A^e{V)/kBT ' ^ ' 

The choice of / as the nonfiuctuating quantity is based both on Fig. [6] and on the fact that 

most experiments are performed under this condition. 

The current has to be sufficiently small in order to yield a voltage proportional to it. 
Our requirement was Jsc ^ O.lekBT/h. This value was inspired by Figs. [5] and [6] and is in 
agreement with experiments 38|, |39|, |40] , that required / < SeksT/h. For experiments that 
used larger currents, the conductivity was current dependent in the interesting range 41 1. 
Our requirement is translated into the condition J < 0.1(1 + a / a')ekBT /h. Since a' was not 
known in advance, we used / = 0.1(1 + a / axhieksT jh, which turned out to be appropriate 
in all cases. In several cases we doubled / and verified that (y) was indeed doubled. 

Figure [7] shows the results obtained for cr', for several values of R and ^/L, while a is 
varied over typically four orders of magnitude. Most of the results sit on a universal curve. 
The line in the graph is 

o'/a = (((Tal/^)-'/' + (4/5) (aAL/a) (40) 

and seems to fit the results reasonably well in this range. The results for ^ = 0.3L are too 
high. We suspect that for this large ^/L ratio the conductivity is infiuenced by the boundary 
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FIG. 7: Superconducting contribution to the electric conductivity, evaluated as described in 
Sec. IIVBI and displayed as a function of the contribution that would be obtained according to 
Eq. (j37p . The line is an empiric function. In all cases A'av = 2 x 10^ and A'j.ciax = 6 x lO"^. For 
CTAL < lOo-, f = 0.003a/-R; otherwise, f = O.OM/R. A R = 5000, ^ = 0.2L, N = 14, TVcont = 2; 
□ ^ = 20, C = 0.2L, N = 14, iVcont = 2; o ^ = 100, ^ = O.IL, N = 28, iVcont = 4; x ^ = 5000, 
^ = 0.3L, N = 14, iVcont = 3. 

conditions; it is known that when ^ is of the order of the length of the wire, qualitatively 
different behavior is obtained 42 1. 



Far from the critical temperature Eq. ( l40l) coincides with the Aslamazov-Larkin result, 
but close to the critical temperature the conductivity increases at a slower rate, even in the 
case P = 0. 



C. Inclusion of the nonlinear term 



Until now we have neglected the term with higher order in ipk- This can be done for 
temperatures that are sufficiently above the critical temperature. We are now interested in 
the region close or below the critical temperature and we therefore set /3 7^ 0. Except for the 
general treatment in Sec. IIV C 11 this section will consider a wire with uniform cross section 
Wq. In the general case, Wq stands for a typical cross section. 
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1. Normalized evolution equations 



We start by defining some useful quantities. and are defined by requiring Pip^wo^js = 
ksT and = h'^/2m. The implication of these requirements is that will be an estimate 

of the size of the order parameter at the critical temperature, and C,/3 will be an estimate of 
the coherence length at this temperature. Solving these two equations gives 

= {2mklT^ / (3hW^f/'' , (41) 
= ih^Wo/APm^kBTY^' . (42) 

In terms of more directly measurable quantities we can write 

= {wo<^l/32n^K^kBTY/' , (43) 

where k is the Ginzburg-Landau parameter; for a wire of cross section wq ~ 10^^^ cm^, 
critical temperature T ~ 1 K and k, ~ 0.1, ■C/s ~ 7 x 10~^ cm. We also define the normalized 
order parameter ip^k = i^k/i^ (and similarly, ip^k = ■ipk/i'), the normalized step size = 
rPip"^ /'yh and the ratio a' = a/l3ip'^ = a(^/3L/L)^; if a is a perfectly linear function of the 
temperature, then a' = (^/3/^(0))^e. 

With these notations, the "macroscopic" increments of ipisk and during a time step r 

are 

^m.ci'/sk = -ffslia' + \^(3k\^)i'f3k + {Nip/Lf{2wk)-^ 

[{wk + Wk+iiii^iik - i^pk+i) + {wk + Wk-i){'ipi3k - i'fsk-i)]] , (44) 

Amacifc = -iff3R/2Wk)i^(3L/Lf[iwo/N){Ih/ekBT) + 

(^/3/L)Im[(wfc_i + Wk)ijj*fjk-ii'i3k + {wk + Wk+iWfiki^pk+i]] ■ (45) 

The numerical details for the integration of these equations are discussed in Appendix |Al 

The fluctuations in the increment of Ak have a variance 2Rfp{wo/Nwk){C,f^l-i/L)'^; the 
fluctuations of ipjBk are taken into account by adding to the real and to the imaginary part 
terms with variance Nfj^WQ^js/wkL and then modifying the phase according to Eq. ( l23l) . 

2. Basic results for low currents 

Figure [8] shows the average value of |'?/'/3fcP, obtained using Eqs. (H^ - (H5|) and Eq. (lASp for 
a wire of uniform cross section and in a region that includes the critical temperature, as a 
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FIG. 8: Average amount of superconducting pairs as a function of q, close to the crit- 
ical temperature. Empty symbols stand for (|'0^p)(a') (lower abscissa) and are fitted by 
Eq. ([Ml); filled symbols stand for (|'i/'/3p)((-^sc)) (upper abscissa) and are fitted by (/gc)^^^ = 

2.61(^^/L)V2(|^^|2)0.58g^^2^/;,, ^ 21, iVcont = 2, iV^v = 4 X 10^ iVrclax = 1-2 X 10^ = 0.2L, 

L = 10, = 5 X 10-1 A ^ = 0.1; □ ^ = 3; ^ = 100. 

and R were varied. iV'/^fcP was averaged over A'^av iteration steps and also over k in the range 
1 + ^cont < k < N — Ncont- Several cases a was swept downwards and then upwards, in 
order to check for hysteresis; for the parameters we studied, no evidence for hysteresis was 
found. As could be expected from Eq. the values of dipisl"^) are a universal function of 
a'. We also evaluated for different values of (.p/L and L (not shown in Fig. [8]), and 

the results lie on the same curve. 

In the considered range, (IV^/jp) can be fitted by the empiric expression 

il^Ppl^) = -1.44(a72 - log,(l + e"'/2)) - 0.70(«' - log,(l + e"')) - 1.93/(1 + e"') . (46) 

As in 



431], our results have a hardly visible inflection point near a' ~ —2.5; our fit 0461) is not 
sufficiently accurate to reproduce this feature. As a test for linear response, we evaluated 
(|'?/'/3p) in several cases for zero current and for I = Q.leksT /h] the current caused a decrease 
of (l^/'^p) of the order of 0.2%. 

We also evaluated (/|c) for zero total current and various values of a', ^/j/i^, L and R. 
Some of the results are shown in Fig. [HI 

Figure [HI shows the superconducting contribution to the electric conductivity for several 
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FIG. 9: Superconducting contribution to the electric conductivity in a region where the order 
parameter is not small. The symbols were evaluated using Eq. (j39p . The line was evaluated using 
Eq. dH]). Unless stated otherwise, N = 21, iVcont = N/7, iVav = 2 x 10^ iVrciax = 6 x lO"^, = 0.2L 
and L = 100. f/j = 10"^ for ^ < 30 and L < 400; otherwise, = 10"^ / = O.lekBT/h, except 
for the two leftmost squares, for which / = 0.2ekBT/h. □ i? = 0.1; A R = 3; • R = 30;+ R = 2, 
L = 200; * ^ = 3, L = 400; x ^ = 3, = O.IL, N = 42; o a' = -1.146 ((|^/3p) = 1); o a' = -2. 

values of R, L and C,i3, with a' in the range —2 < a' < 2. The results were fitted by a 
modification of Eq. (HOl) . Writing = R^^L^ /L^, we defined the quantity 

a/3AL = (irV/8)(a' + (I^^H + 0.3488 tanh(0.734(|7/>^|2)i-5))-3/2 , (47) 

which is a generalization of ctal and reduces to it in the limit ipp — » 0. (iV^/jp) was evaluated 
using Eq. (H6ll . Finally, the fit used for a' was 

^7^ = ((^/3AL/a)-^/^ + 1.205(a^AL/^)-°-'')-'/' . (48) 




3. Voltage distribution 

In the previous sections we have evaluated the time average of the voltage. In this section 
we investigate how the voltage fluctuates over time steps. We define L' = (1 — 2Ncont/N)L, 
the length of the wire excluding the contact regions, and the dimensionless normalized 
voltage V = '-feV/l3ijj'^. For a single step, the normalized voltage drop over L' is = 
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FIG. 10: Distribution of the voltage drop along the wire, evaluated for single steps. L' = 3.2,^^, 
/ = eksT/h, L = W,R = 1,N = 30, iVcont = 3, = 5x10-"^, iVav = 3x 10^ iVrciax = 6x10'^. The 
solid line is the gaussian distribution {27t{{V^) - (F)^))"^/^ exp[-(F - (l7))V(2((y2) _ {V)^))]. A 
a' = —0.4; • a' = —3; o a' = —16. For visibility, the line and symbols for a' = —0.4 (respectively, 
—3) have been shifted 5 units to the right (left). 

Figure [10] is a set of histograms for the values of V obtained in every step, for several 
values of a'. V has contributions from Amac^fc in Eq- (HHl) and also from the fluctuations of 
A^. Since the fluctuations of Aj. are known to be gaussian, the results we present take into 
account the contribution from Amac^fc only. These histograms have been fitted by gaussian 
distributions. 

From Fig. [10] we might conclude that the voltage distribution is approximately gaussian, 
but such a conclusion would be misleading, since the voltage drops for steps that are close 
in time are correlated. A more informative distribution is that of the "consecutive phase 
change," which we define as follows. For every step we check the sign of the macroscopic 
contribution to the phase increment, J2k=i+Ncont '^^k, and add all the consecutive increments 
until the sign is reversed. Note that in order to decide when to stop adding increments and 
start evaluating the next phase change we monitored only the macroscopic contribution to 
^AAfc, but in order to evaluate these increments we took into account the microscopic 
fluctuations as well. The results are shown in Fig. [TTJ Within statistical accuracy, all 
the curves decay monotonically. Although the sizes of typical voltage drops in a single 
step increase when a' becomes more negative, we see that the probability of large phase 
excursions decreases with \a'\. We also see that the probability of large negative phase 
changes is greater than that of large positive phase changes, and this explains the average 
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FIG. 11: Distribution of the absolute values of "consecutive phase changes." Empty symbols are for 
positive changes of the phase and filled symbols for negative changes. The parameters, including 
symbol identity, are the same as in Fig. [TOl 

voltage drop. 

According to the Langer-Ambegaokar-McCumber-Halperin (LAMH) scenario the 
voltage drop is due to phase slips, in which the phase decreases by 2tt. Accordingly, we 
expected to find a relatively large probability density near 27i in Fig. [TT], but we did not 
obtain support for this view. It might seem that this negative result contradicts positive 
evidences, such as the observations by Lukens and Goodkind 4J] or the simulations in 



but these evidences correspond to a closed ring setup, and not to a fixed current. It has 
been further claimed that for a long wire the existence of phase slips ought to be insensitive 



to the boundary conditions (page 1088 in [18|]), but this is not the situation studied here. 



Phase slips will be reconsidered in Sec. IIVC5I 



4- Scaling with length 

In the previous sections we have tacitly assumed that the voltage drop along a uniform 
wire is proportional to its length. We have taken advantage of this assumption and consid- 
ered short wires, that do not require large computational resources. In this section we check 
this assumption. We have evaluated the voltage and its variance for several values of a' and 
several lengths, while intrinsic quantities such as R/L are kept at a fixed value. The results 
are summarized in Table [ITl The average of the voltage per unit length, and its variance per 
unit length, which should in principle be independent of L', seem to increase slightly with 
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TABLE II: Voltage drop and its variance, per unit length of the wire. L' is the length (excluding the 
contact regions) over which V is evaluated. The parameters used are i?/L = 0.033, / = 0.3ekBT/h, 



= 


7.2h/^/2mkBT 


Ncont = 3, iV = 


6 + 3.6L7e^, 


r/3 = 5 X 10-^ iVav 


= 2 X 10®, 


iVrclax = 6 X 10^. 




WHje/kBT)^p{V)/L' 




- {y?)iL' 


a' 


L' = 2.hip 


V = 10^/3 


v = 3oe/3 


V = 2.5^/3 


L' = ioe/3 


L' = 


-0.5 


8.48 


8.87 


8.97 


1.52 


1.58 


1.75 


-1 


5.89 


6.31 


6.42 


2.11 


2.15 


2.23 


-2 


1.26 


1.31 


1.30 


4.35 


4.48 


4.52 



L'; in the case of {V)/L', this trend seems to saturate within the considered range. For the 
parametrers in Table HTl the deviations from the proportionality {V) oc L' are most likely due 
to the finite sizes of the elements of the computational grid and of the contact regions; on 
the other hand, this proportionality breaks down completely when additional length brings 
about an additional phase-slip center [20]. 



5. High- current phase slips 

The phase slips in the LAHM scenario are assumed to be caused by thermal fluctuations. 
We unsuccessfully looked for them in Sec. IIVC3[ On the other hand, early experiments 



2d . l28l | show steps in the volt age- current characteristic, that have been attributed to phase 
slips. These steps appear at currents that are much larger than the current that characterizes 
the scale for thermal fluctuations, eknT/h. Modern versions of these experiments are the 



subject of recent investigations [32, l33|. A numeric study using the same simple TDGL 



Ml 



model as in the present paper was performed by Kramer and Baratoff 

According to the basic phase-slip idea, there are moments at which ip vanishes at places 
called phase-slip centers. When this happens, the order parameter releases part of its winding 
by changing by 27r the phase difference across this center. This phase change is equivalent 
to a change in the electromagnetic potential A, and gives rise to a voltage pulse. According 
to typical models, phase-slip centers are located at fixed positions. 

Figure [12] shows the ratio between the average voltage and the voltage that would be 
obtained for the same current if the wire were in the normal state. The continuous line was 
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FIG. 12: Voltage-current characteristic for high currents. T^iormal stands for L'l /gwq and (• • • ) is 
the average over A'^av steps. The hue (with the hysteresis loop) is obtained when fluctuations are 
ignored and the symbols take fluctuations into account. The inset shows the minimum value of 
the order parameter as a function of the current, a' = —16, L' = 3.2^^, L = 10, i? = 1, = 30, 

iVcont = 3, = 5 X 10-^ iVav = A^rclax = 6 X 10^. 

obtained when thermal fluctuations were ignored, i.e. rj^^v and 7]^ were set equal to zero. 
In this case we find hysteresis. The symbols were obtained using the actual values of riu,v,A 
required by the fluctuation-dissipation theorem. For the parameters used, the inclusion of 
fluctuations leads to results that are independent of history. 

We denote by iV'Imin (and similarly j-?/'/?!!!!!!!) the minimum value of among the elements 
^cont + 1 ^ k < N — Ncont for a given time step. The inset in Fig. [T^ shows, as expected, that 
the region of fast increase of the average voltage is the same as the region of fast decrease 
of (I^A/jImin)- the absence of fluctuations, the smallest values of {ipkl appear always in the 
middle of the wire, but, when fluctuations are included, the positions of the minima are 
widely spread. 

If phase slips give a considerable contribution to the average voltage, then we expect 
that the voltage will be large for steps that have small iV'Imin, since ideally \ip\mm = 
when a phase slip occurs. It follows that V and iV'Imin ^^^^ ^e negatively correlated and 
(^lV'lmin)/(^)(l''Almin) ^^^^ ^c Icss than 1. Figure [T3] shows this ratio as a function of the 
current and we see that there is indeed a region where (^|V'lmin)/(^)(l'^lmin) < 1- 
expect that this is a region where phase slips occur. The position of the minimum of 
(^iV'Imin)/ (^)(l''Mmm) nearly coincides with the current at the lower vertex in the hysteresis 
loop shown in Fig. [121 but the dip extends to currents beyond the hysteresis region. The inset 
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FIG. 13: The ratio (^|V'lmin)/(^)(lV'lmin) measures the correlation between the voltage and the 
minimum size of the order parameter. This ratio has to be less than 1 if phase slips contribute to 
the voltage. The inset shows the region where correlation is significant. The parameters are the 
same as in Fig. [T21 o Fluctuations properly taken into account; ★ size of fluctuations divided by 2. 



FIG. 14: Distribution of consecutive phase changes for Ih/eksT = 92. The parameters are the 
same as in Fig. [T2j The peak near zero has been chopped in order to make phase slips visible, o 
Phase changes restricted to the region above —Sir; ★ phase changes below — Svr were brought into 
the region between —Svr and — vr by addition of the appropriate integer multiple of 27r. 

in Fig. [13] is a close-up at the dip region. In order to appreciate the influence of fluctuations, 
(^l^lmin)/ (^)(lV^lmm) also evaluated with r]u^v,A set to half their actual values. We see 
that fluctuations extend the dip region in both directions, but mainly towards lower currents. 

Figure [H] shows the distribution of consecutive phase changes, as defined in Sec. IIVC3i 
for Ih/eksT = 92. We expect from Figs. [T2l and [T3l that for this current phase slips will be 
present. Most of the consecutive phase changes are close to zero, and have been chopped 
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in the figure. However, there is also a significant amount of cases for which the consecutive 
phase change is close to —271, as predicted. 

For this large current, the sign of the phase change is not necessarily reversed between 
phase slips. Therefore, the consecutive phase change may involve more than one phase slip 
and the consecutive phase change should ideally equal an integer multiple of — 27r. In order 
to take this possibility into account, we added integer multiples of 27r to phase changes that 
were smaller than — Svr, until we brought them to the region between — Stt and — tt. The 
results obtained in this way are shown by the stars in Fig. [Ml 



6. Comparison with experiments 



It is usually accepted that different models are required for the description of the elec- 
tric conductivity in several temperature ranges in the vicinity of the critical temperature. 
Above the critical temperature, conductivity is believed to be dominated by the exchange of 
"cooperons" between normal electrons [45]. Very near the critical temperature, a promising 
method is the Hartree-Fock approximation to TDGL 
but not too near it, the LAMH model gives good fits 
there is a "foot" which in some cases has been attributed to the contacts 
to macroscopic quantum tunneling 31 



21 1 . Below the critical temperature, 
191]. For still lower temperatures 
461 ] and in others 



39| . Moreover, the parameters we have studied in 



section IIVC2I differ by several orders of magnitude from those of the available experimental 
data. In spite of this situation, it is instructive to check how well our naive formalism can 
fit the experimental data. 

We analyzed samples of aluminum 38], indium 39] and tin [4^, in a temperature range 
where lO^^a ^ a' < lO^'^a. We chose the samples so that they could be considered as 
one-dimensional, sufficient data in the analyzed range were available and all the required 
properties of the sample (e.g. length) were reported. We denote by the temperature that 
the experimenters considered to be the critical temperature. The data for aluminum are 
above and the others below T^.. The normal resistance of the tin sample was smaller by 
more than three orders of magnitude than those of the other samples; this was due both to 
its larger cross section and to its longer mean free path. The current passed through the tin 
sample was slightly above the linear range. 



The properties of a sample were estimated as follows. The length was measured by means 
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FIG. 15: Paraconductivity of an aluminum microwire as a function of the temperature above T^. 
The dots are experimental data and the curve is a fit using Eqs. ()46p 



of a microscope; the cross section was either estimated from microscopic observation or from 
the critical current. Assuming the geometry is known, the mean free path can be obtained 
from the normal conductivity The BCS coherence length and the London penetration 
depth are assumed to be the same as in bulk material, in spite of the fact that the critical 
temperature is significantly different than that in bulk. ^(0) and k can be obtained from 



these values and the mean free path [46|], and from Eq. fj43l) . is either assumed equal 
to the critical temperature of a reference sample with large cross section, or is taken from a 
fit to some model. For 7 we take the value Tih? /IQkBT'jn^^^Y . 

We denote by AT = T'^ — the difference between the the critical temperature estimated 
by the experimenters and that required by our model. We fitted the experimental data using 
our phenomenological Eqs. fj46|) - fH8|) . with ^(0), and AT as free parameters. The results 
are shown in Figs. [TSHTGI there is semiquantitative agreement, but the fitting functions are 
too concave. 

The values of the fitting parameters ^(0), and AT are displayed in Table UTTl They 
cannot be taken too seriously, because it is possible to increase or decrease simultaneously the 
three parameters in a wide range without bringing about a clear difference in the goodness of 
the fit; moreover, the values of these parameters depend on the selected range. In all cases, 
Tc is reasonably close to T^,'; for comparison, for the In microwires T^' differs by ~ 800mK 
from the critical temperature in bulk. 

^(0) and enter our equations in two ways. One of them is through a\ and in this case 
only their ratio ^/3/^(0) is of influence. In all cases, the experimental and the fltted ratio 
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FIG. 16: Resistances of indium microwires (upper scale) and a tin whisker (lower scale) as functions 
of temperature. R is the resistance of the sample and is its resistance in the normal state. A 



Sn, tfo ~ 2 X 10 ^cm^; • In, wq ~ 1.7 x 10 -"^^cm^; o In, wq ~ 5.2 x 10 



1 1 2 



-11 



cm 



TABLE III: Parameters used in Eqs. (06]) and 
for the theoretical lines in Figs. [TSHlGi "exp" stands for the values estimated in the experiments 
and "fit" corresponds to the present fit. AT is the difference between the critical temperatures 
estimated in the experiments and those required here. The lengths ^(0) and .^^ are expressed in 
/um. ^(0)°^P is not the same for the thin and the wide In wires, because they have different mean 
free paths. 



Al 
In (thin) 
In (wide) 

Sn 



AT (mK) 
22 
83 
55 
0.355 



^(O)exp 

0.15 
0.042 
0.071 

0.22 



1.1 

0.33 
0.42 

14 



3.1 
1.0 
3.1 

20 



tfit 

5.8 
3 

5.6 

1870 



differ by less than half an order of magnitude, a factor that can be attributed to the coarse 
estimates, both in the experiment and in the fit. ,^(0) or themselves enter through K', 
which in turn depends on the assumption 7 = 7r/),^/16A;BT^'m^(0)^; in this case it seems that 
^(0)*^* and are too big. This might indicate that the theory for 7 is inappropriate and a 
more elaborate theory might be required 481 ]. 
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D. A constriction 



We deal now with the case for which the present formahsm has been designed, namely, 
a wire with nonuniform cross section. Denoting the arclength by s, we consider a wire such 
that its cross section is uniform in segments of length NcontL/N at each extreme, assumes 
the value Wdr 7^ wq at the position s = Sdr, NcontL/N < Sctr < L — NcontL/N, and is 
linear in s for NcontL/N < s < Sdr and Sdr ^ s < L — NcontL/N. We denote by Wq the 
average cross section in the examined region NcontL/N < s < L — NcontL/N, i.e. the cross 
section in the contact regions has to be 2wq — Wctr- If Wctr < Wq there is a constriction, 
and in the opposite case there is a widening. We still use the notations R = ie'^'yL/hawo, 

= {h^WQ/ APm^kBTy/^ . In the discretized version we identify Wi with the cross section at 
s = {i- l/2)L/N. 



1. Conductivity 

With these notations, Eq. (!39|) has to be replaced with 

_ {N - 2Ncont)RMwo/e 2wq - Wctr _ ^ 
a SN^e{V){wo-Wctv) Wctv 

where V is the potential drop excluding the initial and the final Ai'cont elements. 

Figure fT71 compares the conductivities of a uniform and two nonuniform wires, as functions 
of the temperature. As the temperature decreases, the conductivities of the nonuniform 
wires increase at lower rates. An intuitive reason could be that in the thinner parts of the 
wire fluctuations are more effective in breaking superconductivity. Several evaluations of a' 
were repeated for different currents and we verified that our results correspond to the linear 
regime. 



2. Phase slips 

We repeat our study of Sec. IIVC51 but this time the wire has a constriction in the middle. 
The results are shown in Fig. [181 Our results for (V^|V'|mm)/(^) (I'^lmin) indicate that the 
constriction enhances the influence of fluctuations and the phase-slip region is both enlarged 
and shifted to lower currents. A shift to lower currents was also found in Issl. A curious 
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FIG. 17: Contribution of superconductivity to the conductivities of a uniform and two nonuniform 
wires, as functions of the temperature. L = 10, i? = 100, Sctr = 0.5L, = 0.2L, / = O.lefc^T/^, 
N = 21, A^cont = 3, = 5 X 10""^, iVav = 2 X 10^, iVrclax = 6 X 10^. For visibility, the symbols have 
been joined by straight segments. • uniform wire; A Wctr = 0.2?i;o; ^ Wctr = I-Swq. 

effect is that, for currents slightly below this region, the voltage and the minimum size of 
the order parameter are positively correlated. 

The insert shows the probability density for consecutive phase change, for a current in 
the phase-slip region. As expected, the constriction promotes the occurrence of phase slips: 
the probability density near —2tt is much larger than in the case of a uniform wire, and also 
drops faster away from — 27r. Also, when the order parameter vanishes during a phase slip, 
it now happens almost exclusively at the constriction. 

3. Charge accumulation 

As a last test for statistical mechanics, we disconnect the current source and connect a 
capacitor C between the points s = and s = sc- This capacitor will be charged with a 
charge Q that fluctuates in time. The region < s < sc will be forced to sustain a potential 
drop Q/C and the current Ji in this this region will be obtained as a Lagrange multiplier, 
as in Sec. IIII D[ The current in the region sc < s < L has to vanish. During each step, the 
charge in the capacitor increases by the amount Iit. 

From the point of view of the capacitor, the wire is just a heat reservoir. Therefore, 
writing Q = Q/y^ksTC, the probability density for the charge will be proportional to 
exp(-QV2). From here, we expect = 2"/2r((r2 + l)/2)/0F. Table IV] compares this 
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FIG. 18: Normalized voltage and voltage-order parameter correlation as functions of the current. 
■Sctr = 0.5L, WctT = 0.2wq, N = 31; the other parameters as in Fig. [12] (for / < 25ekBT/h we took 
higher values of A^'av). • l^/Vnormai; ^ (^l^lmin)/(^)(l^lmin)- Insert: distribution of consecutive 
phase changes for Ih/eksT = 50, as defined in Sec lIV C"3l o Only includes phase changes restricted 
to the region above —3tt; -k phase changes below —3-k were added after being brought into the region 
between —Sir and — vr, by addition of the appropriate integer multiple of 27r. 

TABLE IV: Average value of according to statistical mechanics and to the method developed 
here. Sctr = 0.36L, sc = 0.48L, C = O-Olj"^ e'^ /ksT; the other relevant parameters are as in Fig. 1171 

n 12 3 4 

2"/2r((n+ l)/2)/0F 0.80 1 1.60 3 

our simulation 0.81 1.02 1.65 3.14 



prediction with the results obtained with our method. 
V. SUMMARY 

We have developed a formalism that takes into account the influence of thermal fluctu- 
ations in a one-dimensional superconductor described by the TDGL model. Our method 
incorporates the fluctuation-dissipation theorem in the form required for integration of a 
discretized stochastic partial differential equation by means of Euler iterations. We have 
pointed out the relationship between the size of the fluctuations and the volume of the el- 
ements into which the sample is divided; to my knowledge, this relationship has not been 



36 



previously addressed in the literature. We have verified that our method complies with 
statistical mechanics. 

Using a gauged order parameter (ip) rather than the raw order parameter (ip), we ob- 
tain a formalism in which the electromagnetic potential is not needed in the macroscopic 
evolution. Special care is required in the use of this transformation; although the thermody- 
namic potential is more naturally expressed in terms of ip, the variable that permits direct 
application of the fiuctuation-dissipation theorem is ip. 

As an application, we have evaluated the contribution of incipient superconductivity 
to the electrical conductivity. Significantly above the critical temperature, we recover the 
Aslamazov-Larkin result. We have focused on the regime immediately below the critical 
temperature. We compared our results with those of old experiments in this regime, and the 
agreement is poor. The lack of agreement was expected and does not arise from the way we 
have treated fiuctuations, but rather from TDGL itself; no available theory reproduces the 
experimental data in this regime (although the Hartree approximation has been successful 
in situations in which fiuctuations may be regarded as effectively one-dimensional 49l|). 

Our method enables us to deal with wires with nonuniform cross section. In particular, 
we studied the influence of fluctuations and constrictions on the presence and behavior of 
phase slips. When a constriction is present, phase slips appear at lower currents. 

Our formalism can be easily generalized to more complex situations, such as multiply 
connected topologies where non-integer fluxes are enclosed, time-dependent applied voltages 
or currents, moderate variations in temperature, or different evolution equations that can 
be cast in the form of dissipative equations. Generalization to models that are not purely 
dissipative (e.g. [5(|) is not immediate. 
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APPENDIX A: NUMERIC CONSIDERATIONS 



A well-known difficulty in the solution of partial differential equations by means of Euler 
iterations is that convergence requires quadratically small time steps as the computational 
grid becomes dense. The standard way for overcoming this difficulty is the Crank-Nicolson 
implicit scheme. However, we did not find that this scheme is useful when stochastic noise 
is present. Our general strategy was to perform analytically the most important parts of the 
iteration and then add a simple Euler iteration for the remaining part. Our final algorithm 
was remarkably stable and enabled us to vary practically all the parameters in our model 
over several orders of magnitude. 



1. Treatment of the inter-cell interaction 

For a dense computational grid ^ ^ 1 and the inter-cell interaction becomes the leading 
term in Eq. (HM . We treat this term separately. Regarding i/) = {ipi, . . . ,iPn) as a vector, 
we have to leading term 

f^MV.. (Al) 

where the matrix M has the elements —2aC,'^/'~fh along the diagonal, aC,'^/'~fh when the 
column and the row number differ by 1, and zeros everywhere else. The ffist and the 
last rows of M depend on the boundary conditions; for instance, for periodic boundary 
conditions Mi^ = M^i = a^'^/'jh, whereas for Neumann conditions Mi^ = M^i = and 
Mil = Mnn = -a^^/-fh. I gnoring the other contributions, Eq. (lAlO can be integrated to 
give 

il;{t + r) = exp(Mr)V(t) . (A2) 

M has eigenvalues that are either negative or zero, so that the largest eigenvalue of exp(Mr) 
is 1. 

We can therefore follow the evolution of ipk by means of steps that consist of three stages. 
The ffist stage is described by Eq. (]A2|) and the other two are as described at the end of 
Sec. IIIDI When considering the contribution of the ffist term in Eq. fl20l) . the leading term 
is dropped, since it has already been taken into account. The matrix exp(Mr) is the same 
for every step, so that it has to be calculated just once. 
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2. Treatment of the condensation-energy term 



The evolution equation (jUj) contains two contributions to AmacV'/3fc: an "internal" con- 
tribution —fjsla' + \ipi3k\'^)ipi3k, that will be considered here, and an "external" contribution 
due to the interaction with the neighboring elements k ±1, that was considered above. 

If only the internal contribution is considered, the evolution equation can be integrated 
and gives V^/3fc(t + r) = (a7[(a'+ |V^/3fc(t)n exp(2a'f^) - |^/i^fc(t)P])i/2^^^(i)_ Substituting this 

expression by a computationally less expensive Pade approximant we obtain the algorithm 
that takes into account the internal contribution. 



3. Step size 

Let us first consider the case /3 = 0. For the validity of the algorithm (fT6|) we require 
T <^ ksT /V^{dG/ dukf . We make the estimates dG/duk ~ awQL^'^\il)k\/N and \il)k\ ~ 
■\J N/awoL. On the other hand, the fluctuations of Jsc, that were estimated in Sec. IIV A 11 
have to be small compared to (/sc)- Therefore, f has to be in the range 

400^/N^,L < f < (A4) 

Similarly, we require r -C ksT /V a^OG / dAkY and (/5c) has to be large compared to the 
Johnson noise. From here we obtain 

800a/A^av^ < f < Na/Ri^ . (A5) 

In addition, we require A^reiax^ ^ 1- 

Our analysis for /? 7^ is similar to that leading to Eqs. ( ]A4I) -( |X5|) . but now we use 
the estimate {ipkl ~ 4'- Eq- (IBSp is still qualitatively true, as shown in Fig. [HI With these 
estimates we obtain 

100^^/N,,L < < (L/N^fsf , (A6) 
SOOiyN^^R^lL^ < < L^/NR£,ll} . (A7) 
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APPENDIX B: AVERAGE SUPERCURRENT 



In this appendix we estimate the characteristic size of supercurrents in the absence of 
voltage. Let us consider a wire with uniform cross section with periodic boundary conditions. 
The average of the supercurrent along the wire is 



2t ^ 



Aea^ Lwq 

'sc = 



J]lm[^;^,_i] , (Bl) 

which can also be written as 



^sc = 2. \^n\ sm ^ . (B2) 

n=0 

The ensemble average of Jsc is zero. Its variance is obtained from the square of Eq. (1B2p . 
2 {Aea^^Lwof sr^ . ^ 12,- ,2 • 27rn . 2Tm' 

The expression for the variance is simplified by substraction of the null quantity (/sc)^ and 
noting that for 72 7^ n' (|(^„p|<^„/p) = whereas (|(^„|^) = 2(|(^„p)^. We obtain 

my ^ (l + 4esin2(7m/iV))2 ■ ^ ' 



For ^ 1, the sum can be replaced with an integral; if in addition we assume ^ ^ 1 
and keep only the leading term, Eq. flB4p becomes 



die) = i^lp . (B5) 
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